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ABSTRACT 


A  two  step  hybrid  perturbation-Galerkin  technique  is  applied  to  the  problem  of  deter¬ 
mining  the  resonant  frequencies  of  one-  or  several-degree(s)-of-freedom  nonlinear  systems 
involving  a  parameter.  In  step  one,  the  Lindstedt-Poincare  method  is  used  to  determine 
perturbation  solutions  which  are  formally  valid  about  one  or  more  special  values  of  the 
parameter  (e.g.  for  small  or  large  values  of  the  parameter).  In  step  two,  a  subset  of  the 
perturbation  coordinate  functions  determined  in  step  one  is  used  in  a  Galerkin  type  approx¬ 
imation.  The  technique  is  illustrated  for  several  one-degree-of-freedom  systems,  including 
the  Duffing  and  van  der  Pol  oscillators,  as  well  as  for  the  compound  pendulum.  For  all  of  the 
examples  considered,  it  is  shown  that  the  frequencies  obtained  by  the  hybrid  technique  using 
only  a  few  terms  from  the  perturbation  solutions  are  significantly  more  accurate  than  the 
perturbation  results  on  which  they  are  based,  and  they  compare  very  well  with  frequencies 
obtained  by  purely  numerical  methods. 


'This  research  was  partially  supported  by  the  National  Aeronautics  and  Space  Administration  under 
NASA  Contract  No.  NAS1-18605  while  the  first  author  was  in  residence  at  the  Institute  for  Computer 
Applications  in  Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23665. 
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1.  INTRODUCTION 


In  this  paper  we  present  and  discuss  a  two-step  hybrid  perturbation-Galerkin  technique 
for  the  computation  of  the  resonant  frequencies  of  nonlinear  oscillating  systems  with  a  finite 
number  of  degrees  of  freedom  involving  a  scalar  parameter  t.  In  previous  papers,  we  have 
developed  and  applied  different  versions  of  the  hybrid  technique  to  several  classes  of  two- 
point  boundary-value  problems  for  ordinary  differential  equations  [2,6,7],  to  some  boundary 
value  problems  for  elliptic  partial  differential  equations  [8],  and  to  some  integral  equations 
of  the  first  kind  which  arise  in  slender  body  theory  [5].  For  each  of  these  classes  of  problems, 
the  method  has  yielded  results  which  are  typically  far  more  accurate  than  the  perturbation 
solutions  on  which  they  are  based,  and  often  provide  at  least  reasonable  solutions  even  for 
values  of  the  expansion  parameter  for  which  the  perturbation  solution  is  meaningless. 

The  idea  of  exploiting  perturbation  expansions  in  conjunction  with  Galerkin  or  variational 
techniques  was  introduced  by  Noor  and  Peters  in  1979  [12]  and  developed  by  Noor  and  his 
collaborators  in  a  number  of  papers  (see  e.g.  [11],  as  well  as  many  other  references  cited  in 
[5]).  Noor’s  “reduced  basis  method”  is  a  combination  of  finite  element  or  other  discretization 
techniques,  perturbation  expansions,  and  Galerkin  (or  variational)  techniques. 

In  general  terms,  our  two-step  hybrid  technique  consists  of  computing  a  few  terms  in  the 
perturbation  expansion  of  the  solution  about  one  or  more  values  of  the  parameter  e  and  then 
using  a  subset  of  these  functions,  with  new  amplitudes,  in  a  Galerkin  type  approximation. 
This  manner  of  combining  the  perturbation  and  Galerkin  approaches  (which  we  will  describe 
in  more  detail  below)  seems  to  overcome  some  of  the  drawbacks  associated  with  each  of  the 
methods  when  they  are  applied  by  themselves,  while  preserving  some  of  the  good  features  of 
each  one.  In  particular,  the  perturbation  method  has  at  least  two  significant  drawbacks.  The 
first  is  that,  for  most  practical  problems,  only  a  few  terms  in  a  perturbation  expansion  can 
be  computed  because  of  the  rapidly  increasing  amount  of  “labor”  that  is  required  to  compute 
each  additional  term.  Secondly,  the  expansion  parameter  must  usually  be  restricted  to  values 
which  lie  close  to  the  point  about  which  the  expansion  was  constructed,  in  order  to  obtain 
approximations  of  acceptable  accuracy.  A  drawback  of  the  Galerkin  method  is  the  problem, 
from  a  practical  point  of  view,  of  selecting  a  small  number  of  “good”  basis  functions.  As  we 
shall  demonstrate  below,  the  functions  determined  by  the  perturbation  method  appear  to  be 
very  effective  basis  functions  and  hence  our  method  overcomes  the  main  drawback  associated 
with  the  Galerkin  method.  Also,  the  new  amplitudes  determined  by  the  method  produce 
approximations  which  are  typically  much  more  accurate  than  the  perturbation  solutions  on 
which  they  are  based. 

In  the  following  sections  we  describe  our  method  in  the  context  of  the  problem  of  de 
termining  the  resonant  frequencies  of  nonlinear  oscillating  systems.  For  simplicity  and  in 
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order  to  illustrate  explicitly  some  of  the  salient  features  of  the  method,  we  consider  first 
systems  with  only  one  degree  of  freedom.  The  method  is  described  in  detail  for  such  sys¬ 
tems  and  then  applied  to  several  model  one-degree-of-freedom  problems.  This  allows  us  to 
obtain  several  explicit  results  and  expressions,  which  provide  insight  into  the  application  of 
the  method  to  more  complicated  systems.  We  then  describe  the  method  in  the  context  of 
a  general  (conservative)  system  with  a  finite  number  of  degrees  of  freedom  which  oscillates 
about  a  stable  equilibrium  state.  The  method  is  then  applied  to  the  classical  problem  of  a 
compound  pendulum.  Observations  about  the  method  are  presented  in  the  final  section. 


2.  ONE-DEGREE-OF-FREEDOM  SYSTEMS 


We  consider  first  the  problem  determining  the  resonant  frequency  v  of  a  one-degree-of- 
freedom  nonlinear  oscillator  which  we  write  in  nondimensional  form  as 


v2  u"  +  u  +  (c/a)  f(a  u,  au  u',  e)  =  0, 

(la) 

u(  0)  =  1,  u'(0)  =  0, 

(16) 

n(x  +  2ff)  =  u(x)  for  all  x  >  0, 

(lc) 

r 2  * 

/  /(a  u,  avu\  t)u  dx  =  0. 

Jo 

(2) 

In  Eq  (la),  a  is  the  maximum  amplitude  of  the  response,  e  is  a  “small”  parameter,  and 
/  is  a  specified  nonlinear  function  of  its  arguments.  (Here  the  dimensional  response  y  of  the 
oscillator  has  been  nondimensionalized  by  defining  u  =  y/a.  Also,  the  primes  in  Eqs  (l)-(2) 
denote  differentiation  with  respect  to  the  nondimensional  time  x  =  u ;t,  where  t  is  time, 
u  =  2n/T  is  the  (unknown)  resonant  frequency  of  the  oscillator,  and  T  is  the  (unknown) 
period  of  the  oscillation.  Then  v  —  w/tco,  where  u>o  is  the  linear  natural  frequency  of  the 
oscillator.)  Conditions  (lb)  and  (lc)  insure  that  u  has  its  maximum  amplitude  at  x  —  0  and 
is  27r  periodic.  Condition  (2)  is  a  necessary  and  sufficient  condition  for  Eq  (1)  to  have  a  27r 
periodic  solution.  We  note  that  in  the  special  case  when  /  does  not  depend  explicitly  on 
u',  i.e.  /  =  /(au,  e),  then  Eq  (la)  represents  a  conservative  system.  Then  condition  (2)  is 
satisfied  automatically  and  Eqs  (1)  can  be  used  to  express  n  as 

v  =  n[j  [F(u,  e)]~l/2  du  , 


(3) 

2  c  rl 

F(u,e)  =  1  —  u2  -f  —  /  /(flj,  c)  dx , 

a  Ju 

where  c  is  defined  by  the  requirements  that  c  >  0  and  F(— c,  c)  =  0.  Thus,  for  this  special 
case,  the  computation  of  n  reduces  to  a  quadrature,  which  will  be  useful  for  comparison 
purposes  with  our  perturbation  and  hybrid  results  (see  [9]). 
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3.  THE  HYBRID  PERTURB ATION-GALERKIN  TECHNIQUE 


We  now  apply  a  slightly  modified  version  of  a  two  step  hybrid  perturbation-Galerkin 
technique  which  has  been  introduced  and  discussed  in  a  series  of  papers  [2,5-8].  In  step  one 
of  the  method,  we  use  the  Lindstedt-Poincare  method  [10]  to  find  approximate  solutions  for 
u,  v,  and  a  in  the  form 

u  =  +  0(eN)t 

j= o 

V  =  1  +  X]  V3^  +  (4) 

J=1 

a  -  ai  e3  +  0(e*)> 

3= 0 

which  are  formally  valid  as  t  — *  0.  In  (4),  the  {uj(x)}  are  the  unknown  perturbation  co¬ 
ordinate  functions,  while  the  constants  {i/j}  and  {a,},  which  determine  the  perturbation 
approximations  to  the  frequency  and  initial  amplitude,  respectively,  are  also  unknown.  To 
determine  these  quantities,  we  substitute  the  expansions  (4)  into  Eqs  (1),  expand  the  result¬ 
ing  equations  in  power  series  in  e  about  e  =  0,  and  then  equate  the  coefficients  of  like  powers 
of  t  on  each  side  of  these  expressions.  In  this  manner,  we  obtain  a  sequence  of  problems  to 
solve  for  each  of  the  unknowns.  In  particular,  we  find  easily  that 

ito(z)  =  cos(x),  (5) 

while  each  of  the  functions  Uj(x)  with  j  >  1  satisfies  a  problem  of  the  form 

u"  +  Uj  =  2  Vj  cos(x)  +  gj,  (6) 


with  Uj(0)  =  u'( 0)  =  0  and  Uj(x  -f  2  7t)  =  «j(x).  Here  g}  depends  on  Uk,  Vk,  and  a*  with 
k  <  j.  Thus,  in  order  for  u3  to  be  2ir  periodic,  the  right  side  of  (6)  must  be  orthogonal  to 
both  cos(x)  and  sin(x),  since  otherwise  “secular”  terms  proportional  to  x  sin(x)  and  x  cos(x) 
will  appear  in  the  solution  for  uj(x).  These  conditions  yield  the  relations 


1  r2  it  r2  n 

i/j  =  —  - —  /  gj  cos (x)dx,  /  g3  sin (x)dx  =  0, 

2  7T  JO  JO 


(7) 


which  are  two  equations  for  the  two  unknowns  i>3  and  a^-j.  In  particular,  for  j  =  1  Eqs  (7) 
become 

1  /2  7T 

iq  =  — - -  /  f(a o  cos(x),  —a0  sin(x),  0)  cos(x)dx, 

2  7T  Go  Jo 

(8) 

r 

0  =  1  f(a0  cos(x),  —  a0  sin(x),  0)  sin(x)dx. 

Jo 
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The  second  equation  in  (8)  is  an  equation  for  a0,  while  the  first  equation  expresses  v\  as  a 
function  of  a0. 

In  the  examples  which  follow,  we  shall  exhibit  several  explicit  expressions  for  the  vari¬ 
ous  terms  in  the  perturbation  expansions  (4),  as  well  as  comment  on  the  accuracy  of  the 
approximations  computed  using  them. 

In  step  two  of  the  hybrid  method,  we  use  a  subset  of  the  perturbation  coordinate  functions 
{uj}  determined  in  step  one  as  both  trial  and  test  functions  in  a  Galerkin  type  approximation. 
In  particular,  we  seek  new  approximate  solutions  2,  2,  and  2  for  u,  v,  and  a,  respectively, 
with 

N- 1 

2  =  n0(x)  +  ^  <5j  u-j(x).  (9) 

j=i 

In  (9),  the  functions  {u;}  are  the  perturbation  coordinate  functions  determined  by  the 
Lindstedt-Poincare  method  in  step  one,  while  the  {(§_,}  represent  new  “amplitudes”  of  these 
coordinate  functions.  (We  note  that  (9)  satisfies  conditions  (lb)  and  (lc)  for  any  choice  of 
the  {6j}.)  To  determine  the  amplitudes  {<§_,},  as  well  as  the  quantities  2  and  2,  we  substitute 
(9)  into  (la)  and  require  that  the  residual  is  orthogonal  to  each  ut t,  0  <  k  <  N  —  1.  Also, 
we  require  that  (2)  is  satisfied  when  u  is  replaced  by  2.  Thus  we  obtain  the  conditions 


/2jt 

/  [i/2  2"  +  u  +  (c/a)  /(a  2,  a  2  2',  e)l  Uk  dx  =  0,  0  <  k  <  N  —  1, 

Jo 

Jr2  n 

'  f(au,auu\e)udt  =  0. 

o 


(10) 


Equations  (10)  are  a  system  of  N -\- 1  equations  to  determine  the  N  + 1  unknowns  6\, . . . , 

2,  and  2.  Although  these  equations  must,  in  general,  be  solved  numerically,  we  note  that 
the  solution  to  this  system  is  a  point  in  ( N  +  1  )-dimensional  space,  where  N  is  reasonably 
small,  while  the  solution  to  Eq  (1)  is  a  continuous  function.  Also,  for  small  values  of  e,  it  is 
reasonable  to  assume  that  the  unknown  quantities  in  (10)  are  “close  to”  the  corresponding 
values  in  the  perturbation  solution  (4)  (e.g.,  S:  Rs  eJ).  Hence,  beginning  with  small  values 
of  c  and  then  proceeding  to  larger  values  of  c,  good  starting  values  are  available  for  the 
unknown  quantities  in  (10),  which  can  be  used  with  an  iterative  method  of  solution,  such  as 
Newton’s  method. 

If  we  set  N  =  1  in  (9),  we  find  that  u  =  u0(j")  =  cos(j-),  while  Eqs  (10)  reduce  to 


22 


c  f2n 

1  4 - —  /  /(a  cos(.r),  —2  2  sin(j),  e)  cos(x)2.r, 

7r  a  Jo 


r  2  n 

/  f(a  cos(a-), 
Jo 


—  a  v  sin(.r),  e)  sin(x)  dx  =  0. 


(ID 
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Equations  (11)  are  a  system  of  two  nonlinear  equations  for  the  two  unknowns  v  and  a. 
In  particular,  if  we  examine  the  solution  to  these  equations  for  small  values  of  e  and  let 
u  =  1  +  V\  t  +  0(e2)  and  a  =  a0  +  0(c),  we  find  from  (11)  that  the  equations  for  iq  and  a0 
are  identical  to  Eqs  (8)  for  tq  and  a0-  Thus,  for  small  values  of  e,  our  hybrid  solution  (with 
N  =  1)  reproduces  the  perturbation  approximations  for  v  and  a,  at  least  to  within  terms 
which  are  0(e2)  and  0(e),  respectively. 

In  the  special  case  when  /  is  independent  of  u' ,  we  see  that  the  second  equation  in  (11) 
is  satisfied  for  any  choice  of  a  and  hence  a  is  an  arbitrary  parameter.  Then  the  first  equation 
in  (11)  yields  the  explicit  expression 

V  =  1  +  —  /  f(a  cos(x),  e)  cos (x)dx  .  (12) 

CL  JO 

In  the  next  section  we  shall  comment  on  the  precise  form  of  Eqs  (10)  for  several  examples, 
as  well  as  comment  on  the  accuracy  of  the  hybrid  approximations  generated  by  their  solution. 

4.  EXAMPLES  -  ONE-DEGREE-OF-FREEDOM  SYSTEMS 

We  now  illustrate  some  of  the  basic  features  of  the  hybrid  solutions  outlined  above  with 
examples  of  several  one- degree- of- freedom  systems.  The  small  parameter  e  which  appears 
in  the  nondimensional  version  of  the  oscillator  Eq  (la)  can  have  several  different  physical 
interpretations.  The  examples  which  follow  illustrate  some  of  these  different  interpretations 
and  hence  indicate  some  of  the  classes  of  problems  for  which  we  feel  the  hybrid  method 
will  be  especially  useful.  We  will  discuss  these  classes  of  problems,  as  well  as  some  possible 
variations  of  the  basic  method,  more  fully  in  the  discussion  section  at  the  end  of  the  paper. 

Duffing  Oscillator 

As  our  first  example,  we  consider  the  Duffing  oscillator  which,  in  dimensional  form,  can 
be  written  as 

V  +u%y  +  ay3  =  0,  y(0)  =  a,  y(0)  =  0,  (13) 

where  u>0  is  the  (linearized)  natural  frequency  of  the  oscillator,  a  is  a  specified  parameter, 
and  the  dots  denote  differentiation  with  respect  to  time  t.  We  let,  T  be  the  (unknown)  period 
of  the  oscillation  and  then  define  u ;  =  2  7r/T,  v  —  u>/u>0  ,  x  =  wt,  and  u  =  yja.  With  these 
definitions,  Eq  (13)  can  be  written  in  the  form  of  Eq  (la)  with  /  =  u3,  i.e. 

v2  u"  +  u  +  e  u3  =  0,  u(0)  =  1,  i/(0)  =  0,  e  =  a  a2/^.  (14) 

Thus,  in  this  example,  e  can  be  interpreted  as  a  measure  of  the  amount  of  nonlinearity  in 
the  restoring  force  in  the  system. 
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For  this  case  condition  (2)  is  satisfied  automatically.  Thus,  the  initial  amplitude  a  is 
arbitrary  and  hence  can  be  incorporated  into  our  expansion  parameter  e.  The  Lindstedt- 
Poincare  method  yields 

Uo  =  cos(x), 

u\  =  —  [cos(x)  —  cos(3x)]  / 32, 

U2  =  [23  cos(x)  —  24  cos(3x)  +  cos(5x)]  /1024, 

(15) 


3  21  ,  81  ,  6549  . 

v  —  1  +  -  e  —  - e2  +  - e3  —  - e4 

8  256  2048  262144 

37737  5  936183  6  7 

+  2097152  6  67108864  C  +  ^  ’’ 

The  hybrid  approximation  Z>  is  obtained  from  Eqs  (10)  which,  in  general,  will  be  cubically 
nonlinear  in  the  amplitudes  {67}.  In  particular,  for  TV  =  1,  Eq  (11)  yields  the  solution 


>+?r 

4  J 


while  for  TV  —  2,  Eqs  (10)  yield  the  relations 


(16) 


v  =  [1  +3e/4-3e^/128  +  3e<5,2/2048]1/2, 
bx  -  e  +  3  e 6t/4  -  9  c b\f 512  +  23  e  <5^/16384  =  0. 


(17) 


The  first  equation  in  (17)  expresses  v  as  a  function  of  S\  and  e,  while  the  second  equation 
is  a  (cubic)  equation  for  S\  as  a  function  of  e.  For  small  values  of  e,  we  see  that  (16)  agrees 
with  (15)  to  within  terms  which  are  0(e2),  while  the  solution  for  v  from  (17)  agrees  with 
(15)  to  within  terms  which  are  0(e3). 

In  Figs  1  and  2,  we  have  plotted  several  approximations  to  u  as  a  function  of  the  parameter 
e.  In  particular,  we  have  plotted  perturbation  approximations  (denoted  by  P[TV]),  where 
TV,  the  number  of  terms  in  expansion  (15),  ranges  from  1  to  6;  hybrid  approximations  ? 
(denoted  by  //[TV])  determined  from  (16)  and  (17)  with  TV  =  1  and  TV  =  2,  respectively; 
and  numerical  solutions  obtained  by  evaluating  the  integral  representation  (3)  for  v.  As  the 
figures  illustrate,  the  perturbation  results  (15)  are  useless  for  values  of  t  above  1  (the  radius 
of  convergence  of  the  perturbation  series),  while  the  one  and  two  term  hybrid  results  are  in 
excellent  agreement  with  the  numerical  results,  even  for  values  of  e  as  large  as  500. 
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Simple  Pendulum 

The  equation  of  motion  of  a  simple  pendulum  can  be  expressed  as 

0  +  w2  sin(0)  =  0,  0(0)  =  a,  0(0)  =  0,  (18) 

with  ojq  =  g/ L,  where  g  is  the  acceleration  due  to  gravity  and  L  is  the  length  of  the  pendulum. 
Here  6(t)  is  the  angle  the  pendulum  makes  with  the  vertical.  Using  the  same  definitions  of 
T,  u,  v,  and  x  as  in  the  first  example  and  letting  u  =  6/a ,  we  find  that  (18)  can  be  written 
in  the  form  of  (la)  as 

v2u"  -1-  u  +  (1/e)  [sin(eu)  —  eu]  =  0,  c  =  a.  (19) 

Thus,  in  this  case,  e  can  be  interpreted  simply  as  the  maximum  amplitude  of  the  oscillation. 
(Condition  (2)  is  again  satisfied  automatically  and  hence  the  initial  amplitude  is  arbitrary.) 
Also,  since  the  nonlinear  term  in  (19)  can  be  expanded  for  small  values  of  e  as  a  power  series 
in  e2,  we  see  that  the  series  in  (4)  involve  only  even  powers  of  e.  In  particular,  using  the 
symbolic  manipulation  system  Mathematica  [15],  the  Lindstedt-Poincare  method  yields  the 
results 

u0  =  cos(x), 

u2  =  [cos(x)  —  cos(3x)] /192, 

u4  =  [17  cos(x)  —  20  cos(3x)  +  3  cos(5x)]  / 61440,  /9n\ 


n  1  2  1  4  23  ,  2519  s  ..  10, 

v  —  1  —  —  e2  + - e4  —  - e6  —  - e8  +  0(e10). 

16  3072  737280  1321205760  v  ’ 

This  expansion  converges  for  |e|  <  7 r. 

The  hybrid  approximations  u  and  u  are  obtained  from  Eqs  (10)  with  (e/a)  f  replaced  by 
[sin(eu)  —  eu]/e.  In  particular,  for  N  =  1,  Eq  (11)  yields 

r  1  r2* . 1 1/2 


■  1  ri,r  1 1/2 

=  —  /  sin(e  cos(x))  cos (x)dx 

1  e  7T  Jo 

=  1 -e2/16  +  e4/1536  +  0(e6),  as  e  -»  0, 


which  agrees  with  (20)  to  within  terms  which  are  0(e4).  In  a  similar  manner,  when  the 
solutions  to  Eqs  (10)  with  N  =  2  are  expanded  for  small  values  of  e,  we  find  that  our 
expression  for  v  agrees  with  (20)  to  within  terms  which  are  0(c6),  and  the  solution  for  62 
is  S2  =  e2  +  (4/16  +  0(e6).  Thus,  our  hybrid  solution  u  (with  N  =  2)  agrees  with  the 
perturbation  solution  for  u  to  within  terms  which  are  0(e4). 
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In  Fig  3  we  have  plotted  various  approximations  to  v  as  a  function  of  e.  Included  are 
perturbation  solutions  (20)  (denoted  by  P[iV])  of  order  0{t2N)  where  N  =  1,2,6;  hybrid 
approximations  (10)  (denoted  by  i/[JV])  where  N  =  1,2  based  on  the  first  N  nonzero  func¬ 
tions  {uj};  and  numerical  solutions  obtained  by  numerically  evaluating  Eq  (3)  for  the  case 
of  the  pendulum  equation.  The  figure  illustrates  that  the  perturbation  method  gives  good 
results  for  e  =  9max  up  to  approximately  e  =  3/t/4.  For  e  near  i r,  the  perturbation  results 
converge  very  slowly  and  give  little  hint  that  the  frequency  must  go  to  zero  as  e  — >  7 r.  In  the 
region  near  e  =  7 7r /8  the  H[ 2]  solution  seems  definitely  superior  to  the  P[2]  solution,  but  for 
larger  amplitudes  the  accuracy  of  the  i/[2]  solution  is  not  good,  though  it  does  hint  that  the 
frequency  is  tending  to  zero  at  e  =  it.  The  //[ 3]  solution  (not  shown)  is  little  better  than 
the  H[ 2]  solution. 

Large  amplitude  oscillations  about  an  unstable  equilibrium 

We  now  consider  the  problem  of  determining  approximations  to  the  frequency  of  “large” 
amplitude  oscillations  about  an  unstable  equilibrium  of  a  nonlinear  system.  As  a  model 
problem  in  this  area,  we  consider  the  (non-dimensional)  problem 

v2  u"  —  u  -fi  u3  =  0,  n( 0)  =  a,  t/(0)  =  0,  (22) 


where  u  is  2w  periodic  and  a  is  now  interpreted  as  the  nondimensional  initial  amplitude. 
The  equilibrium  state  it  =  0  is  unstable  for  the  system  described  by  (22),  while  the  states 
u  =  ±1  are  stable.  However,  stable  oscillations  about  u  —  0  exist  for  a 2  >  2. 

To  study  the  stable  oscillations  about  u  =  0,  we  consider  the  modified  problem 


v2u"  +  it  +  e  (it3  —  2  it)  =  0,  it(0)  =  a,  it ' ( 0 )  =  0, 


u(x  +  2  7r)  =  tt(x), 


(23) 


where  e  is  an  unspecified  parameter.  We  now  make  two  observations.  First,  we  note  that  Eq 
(23)  is  of  the  form  of  Eqs  (1)  with  f  /a  =  it3  —  2  it  and  hence  we  can  construct  approximate 
solutions  which  will  be  formally  valid  as  e  — >  0.  Secondly,  we  note  that,  by  setting  e  =  1  in 
(23),  we  recover  our  original  Eq  (22).  Thus,  our  (formal)  procedure  will  be  to  apply  our  two 
step  hybrid  method,  as  outlined  above,  and  then  set  t  =  1  to  obtain  an  approximation  to  the 
solution  to  (22).  In  this  sense,  we  can  interpret  c  in  this  example  as  a  type  of  “homotopy” 
parameter.  In  particular,  as  c  varies  from  0  to  1,  we  can  think  of  our  problem  (23)  as  varying 
from  a  trivial  problem,  corresponding  to  t  =  0,  which  we  can  easily  solve,  to  the  problem  we 
really  want  to  solve,  i.e.  problem  (22),  corresponding  to  t  =  1. 

Following  the  general  method  outlined  above,  we  first  use  the  Lindstedt- Poincare  method 
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on  problem  (23)  to  obtain 


Uo  =  a  cos(x), 

u\  =  —  (a3/ 32)  [cos(x)  —  cos(3x)] , 
u2  =  (a3/1024)  [(23  a2  —  64)  cos(x) 

—  (24  a2  —  64)  cos(3x)  +  a2  cos(5x)], 


(24) 


+ 


(3  a2  — 

1)  e  - 

'81a6 

63  a4 

2048 

256 

6549  a8 

405  - 

'21a4  3  a2  l\  2 
~256  8  2  J  C 


+ 


+ 


+  £  K  +  0(e5). 


V 262144  2048  '512  16  '  8, 

The  hybrid  approximation  u  is  given  by  (9),  where  the  amplitudes  {<5,}  and  v2  are 
determined  by  Eqs  (10)  with  the  term  (e/a)  /  replaced  by  t  (u3  —  2  u).  In  particular,  setting 
N  =  1  and  e  =  1  in  (11)  we  find 


v  =  [3  a2/ 4  —  1]1//2 

=  ( v/3/2)  |a|  +  0(l/|a|),  as  |a|  — >  oo. 
From  the  integral  representation  (3)  of  //,  we  find  for  this  case  that 

V  =  2^  [/  '  +sin2(1))  -  2]“1/2dx] 

=  (l  +  sm2(x)y1/2dx\  1  +  OU/M) 

=  0.8472  |a|  r  0(l/|a|),  as  a  —>  oo. 


(25) 


(26) 


Thus,  from  (25)  we  see  that  t/  =  0.866  \a\  +  0(l/|a|)  as  |a|  — »•  oo,  which  compares  well  with 
the  exact  asymptotic  result  (26). 

In  Figs  4  and  5  we  have  plotted  various  approximations  to  v  as  functions  of  the  amplitude 
a.  In  Fig  4  it  is  seen  from  the  numerical  solution  that  the  frequency  goes  rapidly  to  zero 
as  a  — >  \/2  from  above.  Also  in  Fig  4  the  perturbation  expansions  (24)  (P[Ar]  for  N  = 
1,2, 4, 8, 16)  indicate  that  the  convergence  of  the  perturbation  solutions  takes  place  only  for 
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a  very  limited  range  of  a  values,  perhaps  y/2  <  a  <  1.75.  The  hybrid  solutions  (//[TV]  for 
TV  =  1,2, 3, 4)  show  rapid  convergence  for  a  >  1.7  and  slower  convergence  nearer  a  =  \/2. 
Even  so  in  the  region  just  above  a  —  \/2,  H[ 4]  seems  much  more  accurate  that  P[16].  Figure 
5  shows  that  the  hybrid  solution,  H[ 2]  continues  to  have  good  accuracy  up  to  a  =  50  and 
that  even  H{  1]  is  a  fairly  good  approximation  for  a  =  50. 

Self  excited  oscillations  -  the  van  der  Pol  oscillator 

As  an  example  of  a  nonconservative  system,  we  now  consider  the  limit  cycle  of  the  van 
der  Pol  oscillator,  for  which  f(au,i/au',t)  =  (a2  u2  —  1  )ai/u'  in  Eq  (la),  i.e. 

v2u"  +  u  +  e(a2u2  -  \)uu‘  =  0,  u(0)  =  1,  u'(0)  =  0,  (27) 

u(x  +  2  7r)  =  u(x)  for  all  x  >  0. 

For  this  case,  t  is  a  “tuning”  parameter  and  can  be  interpreted  as  a  measure  of  the  amount 
of  a  special  kind  of  nonlmear  damping  in  the  system. 

Several  terms  in  the  perturbation  expansions  of  u,  v,  and  a  have  been  reported  [1,4].  In 
particular,  we  have 

rto  —  2  cos(;r), 

u\  =  [3  sin(x)  —  sin(3x)]/4, 

u2  =  [12  cos(x)  —  18  cos(3x)  +  5  cos(5x)]/96, 

(28) 


-  ‘-ir2  + 


=  2  +  ~  e2  — 
96 


—  e4  +  — r6  +  0(es) 
3072  884736  v  ; 


552960 


1019689 

55738368000 


e6  +  0(e8). 


The  radius  of  convergence  for  this  expansion  is  known  to  be  approximately  e  =  1  85  (see 
Andersen  and  Geer  [1]). 

In  Fig  6  we  have  plotted  several  approximations  to  v,  including  the  approximations 
P[N],  obtained  by  including  all  of  the  terms  in  the  perturbation  expansion  (28b)  up  to  those 
which  are  0(e'v+1),  the  approximations  H[N],  which  are  the  hybrid  solutions  based  on  the 
first  TV  perturbation  coordinate  functions  { u} } ,  and  the  approximations  of  Zonnefeld  [16], 
obtained  using  purely  numerical  methods.  The  figure  illustrates  that  the  hybrid  results  are 
superior  to  the  perturbation  results  on  which  they  are  based,  and  appear  to  be  converging 
to  the  numerical  approximations  as  TV  increases.  However,  the  overall  quality  of  the  hybrid 
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approximation  is  not  as  good  as  in  the  previous  three  examples.  In  the  next  section  we  shall 
discuss  a  general  technique  to  improve  the  overall  quality  of  the  hybrid  approximations  and 
then  apply  it  specifically  to  the  van  der  Pol  oscillator. 

5.  COMBINING  TWO  OR  MORE  PERTURBATION  EXPA ?  .dIONS 

In  many  practical  applications,  it  may  be  possible  to  construct  perturbation  expansions  of 
the  solution  u  and  the  frequency  v  about  more  than  one  value  of  the  perturbation  parameter 
e.  That  is,  it  may  be  possible  to  expand  both  u  and  v  as  formal  asymptotic  series  of  the 
form 

u  =  +  0(<*Vh(c)), 

3=0 

(29) 

v  =  5>JaK£)  +  °(Q!PVp+ i(e)), 

j= o 

which  will  be  formally  valid  as  e  — >  ep,  for  p  =  1,2, ...  ,P.  Here  {a?(c)}  is  an  appropriate 
asymptotic  sequence  of  gauge  functions  and  each  up  can  be  determined  completely  by  a 
standard  perturbation  method  (e.g.  a  composite  expansion  of  inner  and  outer  expansions). 
For  example,  in  the  Lindstedt-Poincare  method,  we  have  constructed  expansions  of  the  form 
of  (29)  with  tp  =  t\  =  0  and  a](e)  =  eJ.  In  addition,  it  may  also  be  possible  to  construct 
expansions  of  the  form  of  (29)  which  will  be  valid  as  e  — *  e2  =  oo,  where  the  gauge  functions 
(aj(e)}  might  now  include  terms  such  as  e~?,  where  q  may  be  a  positive  integer  or  a  positive 
rational  number,  and  might  also  involve  appropriate  transcendental  functions,  such  as  log(e) 
(see  e.g.,  [13]). 

A  subset  of  all  of  the  perturbation  functions  up  are  now  chosen  as  coordinate  functions 
for  the  hybrid  technique  and  an  approximation  u  for  u  is  sought  in  the  form  of  (9),  where 
each  Uj  is  now  one  of  the  perturbation  coordinate  functions  up.  To  determine  the  unknown 
amplitudes  we  again  apply  the  conditions  (10). 

To  illustrate  the  application  of  our  method  when  more  than  one  perturbation  expansion 
of  the  solution  is  available  to  us,  we  consider  first  the  problem  of  determining  the  resonant 
frequency  of  a  simple  mechanical  system  consisting  of  a  mass  m  restrained  by  two  identical 
springs,  each  having  natural  length  L  and  spring  constant  k,  midway  between  two  parallel 
walls  a  distance  2d  apart,  as  discussed  by  Arnold  and  Case  [3].  Thus,  if  we  let  h  be  the 
maximum  displacement  of  the  mass,  its  displacement  along  the  centerline  between  the  two 
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planes  is  described  by  hu(x),  where 


v2  u"  +  u  +  n  u  [l  -  (1+-2J^'7?]  =  0,  «(0)  =  1,  «'(0)  =  0, 


u(x  +  2  7T )  =  u(x). 

Here  t  =  h/d,  n  =  A/(l  —  A ),  X  =  L/d  <  1,  v  =  2  7 x/T  cj0,  uo  =  [2  k  (1  —  A )/m]1^2,  and 
a:  =  wo  ^  t,  where  f  is  time.  Thus  e  can  again  be  interpreted  as  a  measure  of  the  maximum 
amplitude  of  the  oscillation  of  the  system. 

For  small  values  of  e,  u  and  v  have  expansions  of  the  form  (29),  with  a*  =  e2U  Using 
the  Lindstedt-Poincare  method,  we  find 

u  =  cos(a;)  +  (/i/64)[cos(3x)  —  cos(x)]  e2  +  0(e4), 

(31) 

v  =  1  +  (3///8)e2  +  0(t4),  as  e  — +  0. 

For  large  values  of  e,  u  has  a  (composite)  expansion  in  the  general  form  of  (29)  with  a2  =  e"-7, 
in  which  the  formal  outer  expansion,  valid  for  jeu|  >  1,  must  be  supplemented  by  inner 
expansions  around  x  =  tt/2  and  x  =  3tt/2,  where  u  vanishes.  In  this  case,  it  is  easy  to  show 
that  u  is  still  well  approximated  by  the  first  term  on  the  right  side  of  (31a),  while 

v  =  1  +  n  +  0(e_1),  as  e  — >  oo.  (32) 

To  apply  our  hybrid  method  to  this  problem,  we  look  for  an  approximate  solution  u  for 
u  in  the  form 

N 

U  =  J2uj(x)  ^0=1,  (33) 

7=0 

where  Uq(x)  —  cos(cc)  and  the  remaining  {uj(x)}  can  be  selected  from  either  the  small  or 
large-e  expansions  of  u.  Substituting  (33)  into  (10),  we  see  that  the  orthogonality  conditions 
become 

Fj{v2,6u62,...,6N)  -0,  j  =0,  l,...,Af, 

/■2?r  r  i  i  (34) 

Fi  =  l  [52«"  +  S  +  /'S(l-(1  +  (252)1/2)j“>^- 

Equations  (34)  are  a  system  of  N  +  1  equations  for  the  N  +  1  unknowns  v2  and  6;,  (j  = 
1,2, ... ,  N).  For  the  special  case  when  N  =  0,  we  can  solve  (34)  explicitly  and  find 


v2  =  1  + 


r2ir  \ 

1  -  77— ~2 - i(~  win  cos2(-v)dx- 

3  (1  +  co s2(a:))1'2 


If  we  can  expand  (35)  for  small  values  of  t,  we  recover  the  expansion  for  v  in  (31b),  while  if 
we  expand  it  for  large  values  of  c  we  recover  the  expansion  (32).  In  Fig  7  wo  have  plotted 


u  determined  by  (35)  as  a  function  of  e,  and  have  also  plotted  some  corresponding  values 
determined  by  the  numerical  evaluation  of  Eq  (3).  As  the  figure  indicates,  the  agreement  of 
the  hybrid  results  with  the  frequencies  determined  from  Eq  (3)  is  excellent. 

As  a  second  example,  we  again  consider  the  van  der  Pol  oscillator.  For  large  values  of  e, 
both  u  and  v  have  expansions  of  the  form  of  (29)  (see  e.g.  [13]),  where  the  gauge  functions 
now  include  e-1,  e~7/3,  e~3log(e), _  In  particular, 


_3-k>g(4)£-1  + 

Z  7T 


The  leading  term  in  the  corresponding  perturbation  solution  for  u  involves  several  inner  and 
outer  expansions,  which  can  be  combined  in  a  straightforward  manner  to  form  the  leading 
term  in  a  composite  expansion  of  the  solution.  We  shall  denote  the  leading  term  in  this 
expansion  by  Uoo(x,e). 

We  now  use  u^x^e)  as  one  of  our  coordinate  functions  in  our  hybrid  solution  (9),  which 
we  write  in  the  form 

N-l 

u  =  80u0  +  ^2  6j  Uj  +  (1  -  60)  Uqo,  (37) 

i= l 

where  {uj,  j  =0,1,2,.. .}  represent  the  perturbation  coordinate  functions  determined  from 
the  (regular)  perturbation  expansion  about  e  =  0.  Thus,  our  solution  (37)  combines  N 
terms  from  the  small-e  perturbation  expansion  of  u  with  the  leading  term  in  the  large-e 
perturbation  expansion  of  u.  The  N  +  2  unknowns  So, . . .  ,Sjv- 1,  a  and  u  are  then  determined 
by  substituting  (37)  into  the  N  - f  2  conditions  (10),  with  N  replaced  by  iV  +  1  and  with  u ^ 
formally  replaced  by  u 00.  In  Fig  8  we  have  plotted  two  perturbation  approximations  to  u\  the 
three-term  small-e  approximation  denoted  by  P0[3]  and  the  one-term  large-e  approximation 
(Eq  (36))  denoted  by  /^[l].  We  also  show  three  hybrid  approximations:  i/[3,0],  which  is 
based  on  the  three  small-e  perturbation  functions;  H[ 0, 1],  which  is  based  on  the  one  large-e 
perturbation  function;  and  H[3, 1],  which  combines  the  information  from  the  three  small-e 
perturbation  functions  with  the  one  large-e  perturbation  function.  As  the  figure  illustrates, 
the  hybrid  approximation  H[ 0, 1]  is  a  considerable  improvement  over  P^,  but  it  gives  poor 
results  for  e  <  2.  The  hybrid  approximation  // [3, 0]  is  a  better  approximation  than  P[3],  but 
it  gives  poor  results  for  e  >  3.  Finally,  //[ 3,  1]  gives  good  results  for  very  small  or  very  large 
values  of  e  and  reasonable  (but  not  accurate)  results  for  “intermediate”  (i.e.  neither  “small” 
nor  “large”)  values  of  t  where  neither  perturbation  expansion  is  accurate. 


6.  SYSTEMS  WITH  SEVERAL  DEGREES  OF  FREEDOM 


We  now  wish  to  study  periodic  solutions  Q(i,i )  to  systems  of  nonlinear  equations  de¬ 
scribing  oscillating  systems  with  several  degrees  of  freedom,  which  can  be  expressed  in  the 


form 


(38) 


AQ  +  B  Q  c  f  (Q ,  Q ,  Q ,  e)  —  0. 

Here  e  is  a  small,  dimensionless  parameter,  the  dots  denote  differentiation  with  respect  to 
time  t,  Q  is  an  M-component  vector  of  dependent  variables  (where  M  >  1  is  the  number 
of  degrees  of  freedom  of  the  system),  A  and  B  are  real,  symmetric,  positive  definite,  M  by 
M,  constant  matrices,  and  /  is  an  M-component  vector.  (Here  it  is  convenient  to  require 
the  components  of  A  to  be  dimensionless,  while  the  components  of  B  and  /  are  required  to 
have  the  units  of  t~2.)  We  assume  that,  for  small  values  of  e,  the  term  t  f  can  be  expanded 
in  a  series  of  the  form 

e/  =  e/i  +  e2  f 2  +  .. .  (39) 

where  each  /,  is  independent  of  e.  We  also  assume  that  the  eigenvalues  \:  of  the  generalized 
eigenvalue  problem  associated  with  the  linearized  version  of  Eqs  (38),  i.e. 

[B  —  Aj  A]  Qj  =  0,  j  =  1,2, . . . ,  M,  (40) 

are  distinct  and  simple,  and  that  the  corresponding  eigenvectors  aj  have  been  normalized 
so  that  ( Aclj ,  aj)  =  1  for  each  j.  (Here  (  ,  )  is  the  usual  Euclidean  vector  inner  product.) 
Consequently,  the  set  of  eigenvectors  {<5q,  1  <  j  <  M)  is  complete  and  A-orthonormal. 
Then,  for  e  =  0,  Eqs  (38)  have  periodic  solutions  of  the  form 

Q(t)  -  a  cos(ujp  t  -  0)ap,  =  \XJ2,  (41) 

where  a  and  (3  are  arbitrary  constants  and  1  <  p  <  M. 

Let  T  be  the  (unknown)  period  of  a  periodic  solution  of  (38)  which  reduces  to  the  solution 
(41)  as  t  — ►  0.  Then  we  define 

u>  =  -=jr,  v  =  — ,  x=ut  =  vuvt,  q(x)  =  Q(t).  (42) 

1  u>p 

In  terms  of  these  variables,  Eqs  (38)  become 

v2  u2  A  q"  +  Bq  +  tf(q,  vupq\  v2  u2p  q",  )  =  0, 

q{x  +  2  7T )  =  q{x),  (43) 

where  the  primes  denote  differentiation  with  respect  to  x. 

In  the  first  step  of  our  hybrid  technique,  we  again  use  the  Lindstedt-Poincare  method 
to  construct  perturbation  solutions  which  are  formally  valid  for  small  values  of  e.  Thus,  we 
seek  solutions  to  Eqs  (43)  in  the  form 

q(x,e)  =  q0(x)  +  tq  i(.r)  +  c2  J)  +  ■••* 
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qo(x)  =  a  cos(x  -  /?)  ap, 

fy  =  1  -|~  t  U\  -b  €*  1*2  “b  •  •  -  • 


(44) 


In  Eqs  (44),  each  qj  is  2rr  periodic  in  x,  and  each  Vj  is  a  constant. 

Substituting  (44)  into  (43)  we  find  that  the  terms  which  are  0(1)  on  the  left  side  of  (43) 
sum  to  zero  (from  the  definition  of  ^0),  while  the  terms  which  are  0(tk)  with  k  >  1  yield  the 
relations 

Wp  A  q'k'  +  Bqk  =2  vku\a  cos(x  -  0)  Aap  +  gk,  k  -  1,2, ,  (45) 

where  gk  involves  qj  and  vj  with  j  <  k.  In  particular,  gx  =  —  fx{qo,  ^P9o>^p9o)-  The 
general  solution  for  qk  will  be  the  sum  of  a  particular  solution  to  Eqs  (45)  and  an  arbitrary 
multiple  of  qo,  which  is  .the  solution  to  the  homogeneous  version  of  Eqs  (45).  To  make  our 
solution  unique,  we  impose  the  orthogonality  condition 

[  {Aqk,  qo)  dx  =  0,  for  k  >  1.  (46) 

Jo 

We  suppose  at  this  point  that  k  is  fixed  and  that  qj  and  i/j  are  known  for  j  <  k.  If  we 
now  take  the  inner  product  of  Eqs  (45)  with  cos(x  —  /?)  ap,  integrate  the  resulting  expression 
with  respect  to  x  between  0  and  2rr,  and  use  the  periodicity  condition  for  qk,  we  find  that 
the  left  side  vanishes  and  hence  we  can  express  i/k  as 

1 

^  =  - T~  l  (9k,  aP)  cos(x  -  P)dx.  (47) 

Z  7T  Up  a  Jo 

In  particular,  for  k  =  1  in  (47)  we  obtain 

\  r2  jt  .  . 

ux  =  - - —  [fx(qo,  ujpqo,  u2  ^),  apJ  cos(x  -  0)dx.  (48) 

Z  7r  up  a  Jo  v  ' 

With  vk  determined  by  (47),  we  look  for  a  solution  for  qk  in  the  form 

Jk  [  M 

9k  =  £  cos(j(x  -  0))  Qp  -  (49) 

3=0  [p=l 

where  the  positive  integer  Jk  and  the  constants  7 ktJtP  are  to  be  determined.  To  determine 
these  quantities,  we  first  express  the  right  side  of  Eq  (45)  as  a  trigonometric  polynomial  in 
the  form 

Jk 

2uku2pa  cos(x  -  0)  Aap  +  gk  -  ^  gk%]  cos {j(x  -  /?)),  (50) 

3= 0 

where  the  vectors  {gk,j}  are  independent  of  x  and  the  inner  product  {gktX,  ap)  =  0.  This 
follows  from  the  definition  of  vk  in  (47).  Substituting  (49)  and  (50)  into  (45),  using  the 
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relations  (40),  and  taking  the  inner  product  of  the  resulting  expression  (40)  with  am,  we  find 
that 

(dk,j ,  5m) 


—  o  -o  o  ' 

m  J  p 


0  <  j  <  Jk, 


1  <  m  <  M,  (j  —  l)2  -f  (m  —  p)2  >  0, 

while  the  constant  7/t,i,p  is  zero  from  condition  (46).  In  deriving  Eq  (51),  we  have  assumed 
that  ^  j2  u2,  unless  j  =  1  and  m  =  p. 

In  the  second  step  of  our  hybrid  method,  we  seek  new  approximations  q  and  u  to  q  and 
i/,  respectively,  where 

N-i 

q  =  q0(x)  +  SjWqA1)-  (52) 

j=i 

In  (52),  the  qj(x)  are  the  perturbation  coordinate  functions  which  appear  in  the  expansion 
(44),  while  the  constants  {<!>j}  represent  new  “amplitudes”,  which  must  be  determined.  We 
note  that  q  is  2 tt  periodic  for  any  choice  of  the  {6j},  since  each  qj  is  2n  periodic.  To  determine 
the  amplitudes  {<5j},  as  well  as  V,  we  substitute  (52)  into  (43)  and  demand  that  the  residual 
is  orthogonal  to  each  of  the  perturbation  coordinate  functions  qk  ,  i.e. 

J0  Aq"  +  Bq  +  tf(q,  Puip|',  P2u>“?",  «)].?*)  ***  =  0,  k  =  0, 1, . . . ,  N  -  1. 

(53) 

Equations  (53)  are  a  system  of  N  equations  for  the  N  unknowns  and  v.  In 

particular,  setting  N  =  1  in  (52)  we  find  that  q  =  q0  while  (53)  yields  the  expression 

£  r2  7r  *1/2 

V  =  1  + - (/{%,  vup%,  v2u2  %),  Qp)  cos(jr  -  0)dx  .  (54) 

7T  cu*  a  Jo  v  ' 

We  note  that,  for  small  values  of  e,  from  (54)  we  can  write  v  =  1  +  tv\  +  0(e2),  where  u\ 
is  given  by  (48).  Thus,  for  small  values  of  e,  our  hybrid  approximation  (with  N  =  1)  agrees 
with  our  perturbation  result  up  to  terms  which  are  0(c2). 


7.  EXAMPLE  -  COMPOUND  PENDULUM 

As  an  application  of  the  results  of  the  previous  section,  we  consider  two  pendula,  each 
of  length  L  and  each  with  attached  mass  m.  They  are  attached  in  such  a  way  as  to  form 
a  classical  compound  pendulum.  Letting  0\(t)  and  02{t)  denote  the  angles  the  two  pendula 
make  with  the  vertical,  we  find  that  6\  and  02  satisfy  the  relations 

2  0i  +  02cos(0]  —  02)  +  02  sin(0i  -  02)  +  sin(0i)  =  0 

(55) 

02  cos(0j  —  02)  +  02  —  Of  sin(0i  —  02)  +  sin(02)  =  0. 


16 


In  (55),  Wg  —  g/L ,  where  g  is  the  acceleration  due  to  gravity.  We  now  let  &i  —  eQi  and 
—  to  write  (55)  in  the  form  of  Eq  (38)  with  e  replaced  by  e2,  where 


B  —  Ua 


Q2  [cos(e(<2i  ~  Q2))  -  1]  +  C  Q\  sin(e(Qi  -  Q2)) 

+‘2c_1 [sin(e  Q\ )  -  cQi] 

Qx  [cos(e(<2i  -  Q2))  -  \]  -  t  Q\  sin(e(<5i  -  Q2)) 
+c_1  u2e  [sin(e  Q2)  -  c  Q2] 


2  -(1/2)Q2(Qi  -  Q2)2  +  QUQ 1  -  Q2)  ~  (1/3 )u:]Q\ 

-(l/2)Qi(Q1  -  Q2f  -  Q\{Qx  -  Q2)  -  (l/6)^|Qi 

For  this  problem  we  find 

U3\  =  —  V2  LJ$ ,  U!2  =  V2-(-  \/2  LVg , 


+  0(e4). 


'  1  ’ 

1 

,V2)  ’ 

°‘  =  C2{-^\ 

c,  =  (l/2)(2  -  v/2),  c2  =  (l/2)(2  +  v/2). 

Thus,  for  small  amplitude  oscillations,  the  periodic  mode  corresponding  to  p  =  1  (i.e.  the 
mode  with  frequency  uq)  represents  an  oscillation  for  which,  at  any  instant  of  time,  the 
pendula  are  always  on  the  same  side  of  the  vertical.  We  shall  refer  to  this  mode  as  the 
“symmetric”  mode  of  oscillation.  In  a  similar  manner,  the  periodic  mode  corresponding  to 
p  —  2  (with  frequency  u2)  represents  an  oscillation  for  which  the  pendula  are  on  opposite 
sides  of  the  vertical,  which  we  shall  refer  to  as  the  “asymmetric”  mode  of  oscillation.  Then, 
setting  p  =  1  in  our  formulae  of  the  previous  section  (with  a  =  1  and  3  =  U  in  the  definition  of 
</o),  we  find  that  we  can  express  the  perturbation  solutions  for  the  frequency  u>  of  oscillation 
of  the  symmetric  mode,  as  w'ell  as  the  angular  displacements  0\  and  02,  as 


U)2/u)g  =  (2  —  \/2)  1 3 


v/2  + 


2  (  —  102  +  71v2)  4  (470602  -  332533 v 2) 

C  1G  +  '  14336  + 


6( -2676206382  +  1 892439443 >/2)  8 

12845056  + 


9\(x )  —  t  q\(x)  =  e  cos(x)  + 


(1  —  \/2)  cos(x)  23(  — 19  +  24\/2)  cos(3x) 
32  +  2688 


(233  +  1375v/2)  cos(x)  |  (812575  -  535236\/2)  cos(3x) 
28672  +  802816  + 


3(970497  -  751760v/2)cos(5x) 
9748480 

02(x)  =  e  q2{x)  =  e  \fl  cos(x)  + 


+  0(t7), 


(2  -  y/2)  cos(x)  (576  -  781^)  cos(3x) 
32  +  2688 


(-2750  -  233\/2)  cos(x)  _  (-1148584  +  756783v/2)  cos(3x) 

28672  +  802816  + 


(-4066960  +  3111 099\/2 )  cos(5x)j  7 

9748480  J  +  ^  ’’ 


It  follows  that  the  perturbation  expansions  for  v  and  for  the  (total)  energy  mug  L2  E  for 
the  symmetric  mode  are  given  by 


v  =  1  +  e 


2  (-31  + 20^)  4  (-61839  +  46888y/2) 


28672 


,  (-284017367  +  197610484\/2)  8. 

+  e - —  +  0(t  ), 

25690112  v  ' 


2  ,  4  3  (29  +  8\/2) 


E  =  2t2  +  c 


6 (2738185  -  550288\/2)  8 

+  t  "  401408  + 

If  6\{x)  and  02{x)  are  computed  to  order  0(e2iV_1),  then  the  energy  E  is  conserved  (in  time) 
to  0{e2N). 

The  corresponding  expressions  for  the  asymmetric  mode  of  oscillation  are  the  \/2  con¬ 
jugates  of  these  expressions,  i.e.  everywhere  \/2  appears  in  Fqs  (59)-(62)  it  is  replaced  by 

-V2. 

For  both  the  symmetric  and  asymmetric  modes  u  is  a  monotonic  decreasing  function  of 
E  in  the  interval  0  <  E  <  6,  the  energy  interval  in  which  the  pendulum  exhibits  a  back- 
and-forth  motion.  For  any  given  energy  in  this  interval  the  frequency  of  the  asymmetric 
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mode  is  higher  than  the  frequency  of  the  symmetric  mode.  Parametric  plots  of  i/(c)  vs.  E(t) 
as  computed  to  orders  0(t2),  0(e4),  0(c6),  0(t8),  0(c10),  and  0(t12),  are  shown  in  Fig  9 
for  the  symmetric  case  and  in  Fig  10  for  the  asymmetric  case.  These  curves  are  labelled 
P[l]  through  P[6].  The  hybrid  results  are  useful  well  past  the  radii  of  convergence  for  the 
perturbation  results.  They  appear  to  be  converging  to  the  correct  values  for  the  energy 
interval  shown. 

8.  DISCUSSION  AND  CONCLUSIONS 

Each  of  the  examples  we  have  considered  provides  some  insights  into  the  hybrid  method 
which  we  now  discuss  briefly. 

The  Duffing  oscillator  illustrates  the  fact  that  the  usefulness  of  a  perturbation  expansion 
can  often  be  limited  because  the  expansion  has  a  finite  radius  of  convergence.  In  this  case, 
the  perturbation  series  converges  only  for  |t|  <  1.  since  the  solution  u.  when  regarded  as  a 
function  of  both  x  and  c,  has  a  singularity  at  t  =  —  1.  (This  follows  from  the  fact  that  the 
solution  to  (14)  is  non-oscillatory  for  e  <  —1.)  Thus,  the  perturbation  expansion  fails  to 
converge  over  the  range  1  <  t  <  oo,  which  is  most  of  the  region  of  physical  interest.  The 
hybrid  method  overcomes  this  limitation  and  provides  meaningful  and,  in  this  case,  at  least, 
very  accurate  approximations  to  the  resonant  frequency,  even  for  c  well  beyond  the  radius 
of  convergence  of  the  perturbation  solution  (see  Figs  1  and  2). 

The  perturbation  expansion  of  the  frequency  (and  solution  u )  for  the  simple  pendulum 
converges  over  the  entire  interval  of  interest,  i.e.  —it  <  c  <  tt,  and  hence,  in  principle,  could 
be  used  to  compute  approximations  to  the  resonant  frequency  for  any  value  of  e  in  this 
range.  However,  the  rate  of  convergence  is  very  slow  for  values  of  t  close  to  ±rr  and  a  great 
many  terms  in  the  perturbation  series  would  be  needed  to  provide  a  solution  of  acceptable 
accuracy  (see  Fig  3).  By  contrast,  the  hybrid  approximation,  based  on  only  two  terms  of 
the  perturbation  expansion,  provides  an  approximation  to  the  resonant  frequency  which  is 
essentially  as  accurate  as  the  perturbation  solution  based  on  six  terms  for  e  not  close  to  ±7 r, 
and  is  reasonably  close  to  the  numerically  computed  frequencies,  even  when  <  is  close  to  ±tt. 
Thus,  we  see  that  the  hybrid  method  has  the  effect  of  accelerating  the  convergence  of  an 
otherwise  slowly  converging  series. 

In  the  example  of  oscillations  about  an  unstable  equilibrium,  we  saw  that  t  he  perturbation 
expansions  depend  on  the  parameter  a  (the  maximum  initial  amplitude),  as  well  as  upon 
c.  Thus,  the  convergence  of  the  perturbation  solutions  will  also  depend  upon  the  parameter 
a.  In  particular,  for  c  =  1,  which  corresponds  to  the  case  of  interest  in  this  problem,  we 
conclude  from  Fig  4  that  the  series  solution  lor  the  resonant  frequency  converges  only  for  a 
very  limited  range  of  values  of  </,  this  range  being  from  a  =  s/2  to  approximately  a  %  1.75. 
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This  is  in  contrast  to  the  hybrid  solutions  which  appear  to  converge  to  the  exact  frequencies 
for  all  values  of  a  >  y/2.  However,  the  rate  of  convergence,  even  for  the  hybrid  solutions,  is 
much  slower  for  values  of  a  close  to  (see  Fig  1)  than  for  larger  values  of  a  (see  Fig  5). 

The  van  der  Pol  oscillator  is  an  example  of  a  nonconservative  system,  for  which  both  the 
frequency  v  and  the  initial  amplitude  a  ot  the  motion  appear  as  unknowns  in  the  problem 
formulation.  For  this  case,  we  again  note  (see  big  6)  that  the  hybrid  solutions  appear  to  be 
converging  to  the  “exact”  (numerically  determined)  frequencies  as  the  number  of  terms  in 
the  approximation  is  increased.  This  convergence  appears  to  hold  even  for  values  of  e  greater 
than  the  radius  of  convergence  of  the  perturbation  expansion,  which  is  approximately  1.85. 
However,  the  rate  of  convergence  of  the  hybrid  solutions  for  c  >  1.85  is  not  as  dramatic 
as  in  the  previous  examples.  Although  we  do  not  as  yet  have  a  good  explanation  for  this 
behavior,  we  feel  (intuitively)  that  it  is  due  to  the  rather  complicated  form  of  the  response 
u  as  e  becomes  large  (see  e.g.  [13]).  It  simply  appears  that  this  behavior  cannot  be  well 
approximated  by  the  perturbation  coefficient  functions  which  have  been  included  so  far  in 
our  hybrid  approximation. 

When  perturbation  expansions  of  the  solution  about  more  than  one  value  of  the  parameter 
t  are  available,  it  has  been  our  experience  that  it  is  usually  advantagous  to  use  at  least  one 
term  from  each  of  these  expansions.  This  is  illustrated  in  the  springs-and-mass  example  as 
well  as  for  the  van  der  Pol  oscillator. 

In  the  springs-and-mass  example,  the  exact  form  of  the  first  term  in  the  small-c  per¬ 
turbation  expansion,  i.e.  u0  —  cos(x),  was  used  in  the  hybrid  solution,  along  with  an 
approximation  to  the  exact  form  of  the  first  term  in  the  large-e  perturbation  expansion, 
which  is  again  cos(x).  In  particular,  we  did  not  include  in  our  hybrid  approximation  the 
various  inner  expansions  which  should  be  combined  with  the  outer  expansion  cos(x)  to  form 
a  uniformly  valid  first  approximation  to  the  solution  when  t  is  large.  Nonetheless,  the  hybrid 
results  are  very  accurate  for  all  values  of  the  expansion  parameter  c,  as  well  as  for  an  entire 
range  of  values  of  the  parameter  fi  (see  Fig  7).  This  illustrates  the  fact  that  apparently  one 
does  not  need  to  determine  the  “complete”  form  of  the  perturbation  coefficient  functions 
and  that  merely  a  good  approximation  to  these  functions  may  suffice  for  use  in  the  hybrid 
approximation. 

For  the  van  der  Pol  oscillator,  the  exact  form  of  the  large- 1  perturbation  solution  is  very 
complicated,  although  several  terms  in  the  various  inner  and  outer  expansions  necessary 
to  construct  this  expansion  have  been  computed  [I3j.  However,  in  our  approximation,  we 
used  only  the  leading  term  from  this  expansion  (which  is  “difficult"  to  compute),  along  with 
several  of  the  small-c  perturbation  coordinate  functions  (which  are  "easy”  to  compute)  to 
form  our  hybrid  approximation.  As  illustrated  in  Fig  8,  there  is  a  definite  improvement 
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in  the  quality  of  the  hybrid  approximation,  although  the  agreement  with  the  numerically 
determined  frequencies  is  by  no  means  entirely  satisfactory.  Presumably,  the  accuracy  of  our 
approximations  could  be  improved  by  including  more  terms  from  either  the  small-  or  large-e 
perturbation  expansions.  However,  we  have  not  carried  out  these  calculations. 

The  compound  pendulum  is  an  example  of  the  application  of  our  method  to  systems 
with  more  than  one  degree  of  freedom.  We  note  from  Figs  9  and  10  that  the  perturba¬ 
tion  expansions  of  the  different  resonant  frequencies  appear  to  have  significantly  different 
radii  of  convergence.  In  particular,  the  perturbation  expansion  of  the  higher  resonant  fre¬ 
quency,  corresponding  to  the  asymmetric  mode  of  oscillation,  has  a  much  smaller  radius  of 
convergence  (approximately  t 2  =  0.092)  than  the  expansion  of  the  lower  resonant  frequency, 
corresponding  to  the  symmetric  mode  of  oscillation  (approximately  t2  —  1.2).  Hence,  the  twro 
perturbation  expansions  will  be  useful  for  calculating  approximations  to  the  corresponding 
resonant  frequency  over  different  ranges  of  values  of  energy  E  (approximately  0  <  E  <  0.17 
for  the  asymmetric  case  vs.  0  <  E  <  2.3  for  the  symmetric  case).  In  both  cases,  however, 
the  hybrid  approximations  based  on  only  a  few  perturbation  functions  appear  to  converge  to 
the  numerically  determined  solutions,  even  for  values  of  t  well  beyond  the  respective  radius 
of  convergence.  We  are  currently  developing  the  hybrid  method  in  the  context  of  general 
Hamiltonian  or  Lagrangian  systems  (and  even  more  general  systems  whose  solution  can  be 
characterized  in  terms  of  a  variational  principle)  and  are  applying  the  methods  to  several 
other  examples.  The  results  of  these  investigations  will  be  reported  elsewhere. 

We  note  that  the  perturbation  parameters  t  employed  in  the  above-mentioned  problems 
have  a  variety  of  interpretations,  but  in  each  case  e  =  0  corresponds  to  simple  harmonic 
motion. 

While  we  have  for  reasons  of  space  limitation  only  reported  on  frequency  functions,  the 
hybrid  technique  appears  to  converge  to  correct  mode  shapes  as  well  as  to  correct  frequencies. 

The  hybrid  method  is  more  successful  in  some  of  the  examples  reported  here  than  in 
others.  We  think  it  is  particularly  effective  for  large  values  of  e  in  Fig  2,  for  large  values 
of  a  in  Fig  5,  and  for  all  values  of  n  in  Fig  7.  It  seems  to  be  less  effective,  as  in  the 
simple  pendulum  problem,  when  “directly”  approaching  a  singularity  rather  than  merely 
approaching  a  radius  of  convergence. 

We  also  note  that  the  equations,  such  as  (10)  and  (53),  which  determine  the  hybrid 
amplitudes  {<5j}  and  u  for  nonlinear  systems  will,  in  general,  be  nonlinear  and  have  multiple 
solutions.  Some  of  these  solutions  may  be  ruled  out  on  the  basis  that  they  involve  complex 
numbers.  Of  the  remaining  solutions,  it  appears  that  the  only  ones  of  interest  are  those  for 
which  the  6j(e)  coincide  with  the  gauge  functions  of  the  perturbation  expansion  in  the  limit 
f->0.  Therefore  following  a  solution  path  starting  at  e  =  0  seems  to  be  an  essential  part  of 
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the  method.  In  our  experience  this  leads  in  each  case  to  a  unique  solution. 

It  is  interesting  to  note  that  the  hybrid  method  provides  a  nice  supplement  to  some  of  the 
methods  and  ideas  presented  by  Van  Dyke  [14]  for  improving  the  usefulness  of  a  perturbation 
expansion.  In  particular,  Van  Dyke’s  ideas  are  applicable  when  a  fairly  large  number  of  terms 
in  a  perturbation  expansion  can  be  computed  (usually  with  the  aid  of  numerical  or  symbolic 
computation).  The  coefficients  in  the  series  are  then  analyzed  to  help  uncover  some  of  the 
analytical  structure  of  the  solution  and  then  this  information  is  used  to  recast  the  series 
into  a  different  form  which,  in  general,  will  be  valid  for  a  wider  range  of  parameter  values. 
In  contrast,  the  hybrid  method  requires  only  a  few  terms  (often  only  one  or  two  terms) 
in  the  perturbation  expansion  to  be  computed  and  then  uses  these  terms  to  construct  an 
“improved”  approximation  to  the  solution. 

In  conclusion,  it  appears  that  the  hybrid  perturbation-Galerkin  technique  is  a  useful  way 
to  enhance  the  usefulness  of  perturbation  solutions  to  resonant  frequency  calculation  prob¬ 
lems.  In  particular,  the  hybrid  solutions  v  provide  useful  approximations  to  the  resonant 
frequencies  for  the  examples  illustrated  here,  even  for  parameter  values  well  beyond  the 
radius  of  convergence  of  the  perturbation  solutions  on  which  they  are  based.  We  are  cur¬ 
rently  investigating  application  of  the  method  to  more  general  Hamiltonian  and  Lagrangian 
systems  with  a  finite  number  of  degrees  of  freedom,  as  well  as  to  systems  with  an  infinite 
number  of  degrees  of  freedom,  i.e.  to  partial  differential  equations.  The  initial  results  of  these 
investigations  have  been  very  promising. 

While  there  are  still  many  theoretical  questions  to  be  answered  about  the  hybrid  method, 
it  is  based  on  an  intuitively  plausible  idea,  it  is  relatively  simple  to  implement,  and  it  appears 
to  provide  reasonable  and  often  very  accurate  approximate  solutions. 
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FIG  1.  Comparison  of  hybrid  (solid  lines)  frequency  computations  for  the  Duffing  oscillator 
with  perturbation  (dashed  lines)  and  numerical  (circles)  computations  for  a  small  range  of  e 


FIG  2.  Comparison  of  hybrid  (solid  lines)  frequency  computations  for  the  Duffing  oscillator 
with  numerical  (circles)  computations  for  a  large  range  of  c  values. 


FIG  3.  Comparison  of  hybrid  (solid  lines)  frequency  computations  for  the  simple  pendu¬ 
lum  with  perturbation  (dashed  lines)  and  numerical  (circles)  computations.  The  amplitude 
ranges  from  0  to  n. 


FIG  4.  Comparison  of  hybrid  (solid  lines)  frequency  computations  with  perturbation  (dashed 
lines)  and  numerical  (dots)  computations  for  a  small  range  of  t  values. 
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FIG  5.  Comparison  of  hybrid  (solid  lines)  frequency  computations  with  numerical  (dots) 
computations  for  a  laxge  range  of  e  values. 


FIG  6.  Comparison  of  hybrid  (solid  lines)  frequency  computations  for  the  limit  cycle  of  the 
van  der  Pol  oscillator  with  perturbation  (dashed  lines)  and  numerical  (circles)  computations 
for  a  range  of  e  values. 


FIG  7.  Comparison  of  hybrid  (solid  lines)  frequency  computations  as  a  function  of  e  for  the 
springs-and-mass  system  with  numerical  (circles)  computations  for  a  variety  of  fi  values. 


FIG  8.  Comparison  of  hybrid  (solid  lines)  frequency  computations  as  a  function  of  e  for  the 
limit  cycle  of  the  van  der  Pol  oscillator  with  perturbation  (dashed  lines)  and  numerical  (cir¬ 
cles)  computations.  The  H[ 3,  1]  solution  is  based  on  two  different  perturbation  expansions. 
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FIG  9.  Comparison  of  hybrid  (solid  lines)  frequency  computations  as  a  function  of  E  for  the 
“symmetric”  mode  of  the  double  pendulum  with  perturbation  (dashed  lines)  and  numerical 
(dots)  computations. 


E 

FIG  10.  Comparison  of  hybrid  (solid  lines)  frequency  computations  as  a  function  of  E  for  the 
“asymmetric”  mode  of  the  double  pendulum  with  perturbation  (dashed  lines)  and  numerical 
(dots)  computations. 
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